# =============================================================================================================================
# Sagarzazu/Kl�ver (PSRM): Coalition governments and party competition: Political communication strategies of coalition parties
# Replication Script for Figure 2
# 11.05.2015
# =============================================================================================================================


#install.packages(c('Design','Hmisc','nlme','Mattrix','lattice','car','arm','memisc','mice','Zelig','Amelia','glmmADMB','plotrix'))
library(Design)
library(Hmisc)
library(foreign)
library(nlme)
library(lattice)
library(Matrix)
library(lme4)
library(car)
library(arm)
library(memisc)
library(mice)
library(Zelig)
library(Amelia)
library(glmmADMB)
library(plotrix)      
options(scipen=3)     





# Loading coefficients and variance covariance matrix estimated in STATA
# ======================================================================

rm(list=ls(all=TRUE))
setwd("C:/Users/Heike.Kluever/Dropbox/Kluver&Sagarzazu/PSRM Submission/Replication materials")
fixeff_raw <- read.csv("fixeff.csv")
fixeff_untransposed <- fixeff_raw[1,]
fixeff <- t(fixeff_untransposed)
varcov_raw <- read.csv("varcov.csv")
varcov <- varcov_raw[1:34,]



# 1. Get b and V
# ==============
(b <- as.vector(fixeff))		  
(V <- as.matrix(varcov))		      
(sqrt(diag(V)))                           



# 2. Set up MVN distribution
# ==========================
nsim <- 1000
mvn <- mvrnorm(nsim, mu=b, Sigma=V)          



# 3. Get coefficients
# ===================

dat <- read.dta("Rdata.dta")
attach(dat)
str(dat)



# 4. Computing predicted values
# =============================

yhat <- list(rep(NA,48))
for (i in 1:48){
yhat[[i]] = (
              mvn[,1] * i
            + mvn[,2] * i^2
            + mvn[,3] * mean(ep_election)
            + mvn[,4] * mean(reg_election)
            + mvn[,5] * mean(mip)
            + mvn[,6] * mean(unempl)
            + mvn[,7] * mean(bse)
            + mvn[,8] * mean(cdu_crisis)
            + mvn[,9] * mean(terror)
            + mvn[,10] * mean(afghan)
            + mvn[,11] * mean(flood)
            + mvn[,12] * mean(worldcup)
            + mvn[,13] * mean(crisis)
            + mvn[,14] * mean(lag_dv)
            + mvn[,15] * mean(issue_dum2)
            + mvn[,16] * mean(issue_dum3)
            + mvn[,17] * mean(issue_dum4)
            + mvn[,18] * mean(issue_dum5)
            + mvn[,19] * mean(issue_dum6)
            + mvn[,20] * mean(issue_dum7)
            + mvn[,21] * mean(issue_dum8)
            + mvn[,22] * mean(issue_dum9)
            + mvn[,23] * mean(issue_dum10)
            + mvn[,24] * mean(issue_dum11)
            + mvn[,25] * mean(issue_dum12)
            + mvn[,26] * mean(issue_dum13)
            + mvn[,27] * mean(issue_dum14)
            + mvn[,28] * mean(issue_dum15)
            + mvn[,29] * mean(issue_dum16)
            + mvn[,30] * mean(issue_dum17)
            + mvn[,31] * mean(issue_dum18)
            + mvn[,32] * mean(issue_dum19)
            + mvn[,33] * mean(issue_dum20)
            + mvn[,34] * 1
            )
}

str(yhat)



# 5. Transform list into matrix
# =============================

yhat_matrix <- do.call(cbind, yhat)



# 6. Taking the means and quantiles of simulated values
# =====================================================

yhat_mean <- list(rep(NA,48))
yhat_lb <- list(rep(NA,48))
yhat_ub <- list(rep(NA,48))

for (i in 1:48){
 yhat_mean[[i]] <- mean(yhat_matrix[,i])
 yhat_lb[[i]] <- quantile(yhat_matrix[,i], c(.050))
 yhat_ub[[i]] <- quantile(yhat_matrix[,i], c(.950))
}

mean <- exp(unlist(yhat_mean))
lb <- exp(unlist(yhat_lb))
ub <- exp(unlist(yhat_ub))
xval <- 1:48

results <- data.frame(xval,mean, lb, ub)



# 8. Plotting predicted probabilites with 90% CI in reversed order of months
# ==========================================================================

xval_rev <- 48:1
orig_data <- data.frame(xval_rev, ub, lb)

dd <- transform(orig_data)

rev_data <- dd[ do.call(order, dd) ,]
rev_data_ext <- data.frame(rev_data, 48:1)


pdf(file="figure2.pdf", height=8, width=11)
par(mfrow=c(1,1),mar=c(5,5,1.5,3))
plot(xval, mean, type="l", col="black", lwd=2, xlab="Months until next election",
      ylab="Predicted diversity in issue attention", xlim=c(50,0), cex=1.5,ylim=c(.015, .04),cex.axis=2.0, cex.lab=2.5)
      lines(rev_data_ext$X48.1, rev_data_ext$ub, col="black", lty=2, lwd=2)
      lines(rev_data_ext$X48.1, rev_data_ext$lb, col="black", lty=2, lwd=2)
      legend(20, .0375, bty = "n", "Point estimate", lwd=2, col="black", cex=1.5)
      legend(20, .0365, bty = "n", "90% Confidence interval", lwd=2, lty=2, col="black", cex=1.5)
dev.off()

